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Abstract 

The decay pattern of aftershocks in the so-called 'coherent-noise' models [M. E. J. 
Newman and K. Sneppen, Phys. Rev. E54, 6226 (1996)] is studied in detail. Ana- 
lytical and numerical results show that the probability to find a large event at time t 
after an initial major event decreases as t~ T for small t, with the exponent r ranging 
from to values well above 1. This is in contrast to Sneppen und Newman, who 
stated that the exponent is about 1, independent of the microscopic details of the 
simulation. Numerical simulations of an extended model [C. Wilke, T. Martinetz, 
Phys. Rev. E56, 7128 (1997)] show that the power-law is only a generic feature of 
the original dynamics and does not necessarily appear in a more general context. 
Finally, the implications of the results to the modeling of earthquakes are discussed. 



1 Introduction 



Dynamical systems which display scale-free behaviour have attracted much 
interest in recent years. Equilibrium thermodynamic systems do only exhibit 
scale-free behaviour for a subset of the parameter space of measure zero (the 
critical values of the parameters). Nevertheless, in nature scale-free systems 
can be found in abundant variety (earthquakes [1], avalanches in rice-piles [2], 
infected people in epidemics [3], jams in Internet traffic [4], extinction events in 
Biology [5] , life-times of species [6] and many more. See also [7] and the refer- 
ences therein) . The origin of this abundance lies probably in the broad variety 
of systems far from equilibrium that can be found in nature. With the onset 
of non-equilibrium dynamics, new mechanisms come into play which seem to 
make scale-free behaviour a generic feature of many systems. However, unlike 
equilibrium thermodynamics, where scaling is thoroughly understood [8,9], for 
non-equilibrium dynamical systems there does not yet exist a unified theory of 
scale- free phenomena (apart from non-equilibrium phase transitions). There 
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do, however, exist several distinct classes of systems with generic scale-free 
dynamic. 

One of the first ideas to explain scale-free behaviour in a large class of dynam- 
ical systems was the notion of Self- Organized Criticality (SOC) proposed by 
Bak, Tang and Wiesenfeld in 1987 [10,11]. They proposed that certain systems 
with local interactions can, under the influence of a small, local driving force, 
self-organize into a state with diverging correlation length and therefore scale- 
free behaviour. This state is similar to the ordinary critical state that arises 
at the critical point in phase transitions, although no fine-tuning of parame- 
ters is necessary to reach it. Since 1987 literally thousands of research papers 
have been written concerning SOC (for a review see [12]), and many different 
dynamical systems have been called SOC (e.g. [13-17]), mostly just because 
they showed some power-law distributed activity patterns. Recently [18] it has 
become clear that several SOC models (sandpile models, forest-fire models) 
can be understood in terms of ordinary nonequilibrium critical phenomena 
(like e.g. [19]). Driving rate and dissipation act as critical parameters. The 
critical value, however, is 0. Therefore, it suffices to choose a small driving 
rate and dissipation to fine-tune the system to the critical point, and this 
choice is usually implicit in the definition of the model. 

Scale-free behaviour does not, however, depend crucially on some sort of criti- 
cal phenomenon. A simple multiplicative stochastic process (MSP) of the form 

x(t + 1) = a(t)x(t) + b{t) , (1) 

where a(t) and b(t) are random variables, can produce a random variable 
x(t) with a probability-density function (pdf) with power-law tail [21-25]. In 
processes of this type, the power-law appears under relatively weak condi- 
tions on the pdf's of a(t) and b(t), thus making the intermittend behaviour a 
generic feature of such models. Applications of Eq. (1) can be found in pop- 
ulation dynamics with external sources [25], epidemics [3], price volatility in 
economics [26], and others. 

Another class of models with a very simple and robust mechanism to pro- 
duce scale-free behaviour has been introduced recently by Newman and Snep- 
pen [27]. These so called 'coherent-noise' models consist of a large array of 
'agents' which are forced to reorganize under externally imposed stress. In 
their simplest form, coherent-noise models do not have any interactions be- 
tween the agents they consist of, and hence, certainly do not display criticallity. 
Nevertheless, these models show a power-law distribution of the reorganization 
events with a wide range of different exponents [28], depending on the spe- 
cial implementation of the basic mechanism. Moreover they display power-law 
distributions in several other quantities, e.g., the life-time distribution of the 
agents. These models have been used to study earthquakes [27], rice piles [27], 
and biological extinction [29-32] . 
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Coherent-noise models have a feature that usually is not present in SOC mod- 
els and is never present in MSP's, which is the existence of aftershocks. In 
most coherent-noise models the probability for a big event to occur is very 
much increased immediately after a previous big event and then decays with 
a power-law. This leads to a fractal pattern of events that are followed by 
smaller ones which themselves are followed by even smaller ones and so on. In 
most SOC models and all MSP's, on the contrary, the state of the system is 
statistically identical before and after a big event. Therefore in these models 
no aftershocks are visible. 

The existence or non-existence of aftershock events should be easily testable 
in natural systems. This could provide a means to decide what mechanism is 
most likely to be the cause for scale-free behaviour in different situations [28]. 
But to achieve this it is important to have a deep understanding of the decay- 
pattern of the aftershock events. 

The goal of the present paper is to investigate in detail the aftershock dynam- 
ics of coherent-noise models. We concentrate mainly on the original model 
introduced by Newman and Sneppen because there can be obtained several 
analytical results. We find a power-law decrease in time of the aftershocks' 
probability to appear, as has been found already in [28]. But unlike stated 
there, we can show that the exponent does indeed depend on the microscopic 
details of the simulation. We find a wide range of exponents, from to values 
well above 1, whereas in [28] the authors report only the value 1. 



2 The model 



We will now describe the model introduced by Newman and Sneppen [27]. 
The system consists of N units or 'agents'. Every agent i has a threshold Xi 
again external stress. The thresholds are initially chosen at random from some 
probability distribution Pthrcsht^)- The dynamics of the system is as follows: 

(1) A stress i] is drawn from some distribution Stress (jl)- All agents with 
Xi < i] are given new random thresholds, again from the distribution 

PthreshV^J- 

(2) A small fraction / of the agents is selected at random and also given new 
thresholds. 

(3) The next time-step begins with (i). 

Step (ii) is necessary to prevent the model from grinding to a halt. Without 
this random reorganization the thresholds of the agents would after some time 
be well above the mean of the stress distribution and average stress could not 
hit any agents anymore. 
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The most common choices for the threshold and stress distributions are a 
uniform threshold distribution and some stress distribution that is falling off 
quickly, like the exponential or the gaussian distribution. Under these con- 
ditions (with reasonably small /) it is guaranteed that the distribution of 
reorganization events that arises through the dynamics of the system will be 
a power-law. 

There are several possibilities to extend the model to make it more general, 
without loss of the basic features. Two extensions that have been studied by 
Sneppen and Newman [28] are 

• a lattice version where the agents are put on a lattice and with every agent 
hit by stress its nearest neighbours undergo reorganization, even if their 
threshold is above the current stress level. 

• a 'multi-trait' version where, instead of a single stress, there are M different 
types of stress, i.e. the stress becomes a M-dimensional vector r\. Accord- 
ingly, every agent has a vector of thresholds Xi. An agent has to move in 
this model whenever at least one of the components of the threshold vector 
is exceeded by the corresponding component of the stress vector. 

An extension that is especially important for the application of coherent noise 
models to biological evolution and the dynamics of mass extinctions has been 
studied recently by Wilke and Martinetz [32]. In biology it is not a good 
assumption to keep the number of agents (in this case species) constant in 
time. Rather, species which go extinct should be removed, and there should 
be a steady regrowth of new species. In a generalized model, the system size 
is allowed to vary. Agents that are hit by stress are removed from the system 
completely, but at the end of every time-step a number AN of new agents is 
introduced into the system. Here AiV is a function of the actual system size N, 
the maximal system size iV max and some growth rate g. Wilke and Martinetz 
have studied in detail the function 

NN e 9 

AN = W --N, (2) 

A max + N(e 9 - 1) 1 ; 

which resembles logistic growth. In the limit g — > oo their model reduces to 
the original one by Newman and Sneppen. In the following we will refer to the 
original model as the 'infinite-growth version' and to the model introduced by 
Wilke and Martinetz as the 'finite-growth version'. 



3 Analysis of the aftershock structure 



We base our analysis of the aftershock structure on the meassurement-procedure 
proposed by Sneppen and Newman [28]. They drew a histogram of all the times 
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whenever an event of size > some constant s± happened after an initial event 
of size > some constant So, for all events > Sq. Consequently, we measure the 
frequency of events larger than s\ occuring exactly t time-steps after an initial 
event larger than s , for all times t. This means that we consider sequences of 
events in the aftermath of initial large events. Normalized appropriately, our 
measurement gives just the probability to find an event > s\ at time t after 
some arbitrarily chosen event > sq. For this to make sense in the context of 
aftershocks we usually have s Q > s x . Throughout the rest of this paper we use 
s and s 1 as percentage of the maximal system size. Therefore a value s = 1 
for example means that we are looking for initial events which span the whole 
system. 

We will denote the probability to find an event of size s > s± at time t after 
a previous large event by P t (s > Si). In order to keep the notation simple 
we omit the constant s . It will be clear from the context what s we use in 
different situations. Note that P t (s > si) is not a probability distribution, but 
a function of time t. Therefore, every increase or decrease of P t (s > Si) in 
time will indicate correlations between the initial event and the subsequent 
aftershocks. For t — > oo we expect all correlations to disappear, and hence 
Pt(s > si) to tend towards a constant. 

It is possible to obtain some analytical results about the probability Pt(s > si) 
if we restrict ourself to the model with infinite growth and a special choice for 
the threshold and stress distributions. If not indicated otherwise, throughout 
the rest of this section we assume pthresh(^) to be uniform between and 1, 
and the stress distribution to be exponential: p s tress{v) — exp(— rj/a)/a. 

Furthermore, we focus on the case So = 1. That means that we are looking 
at the events in the aftermath of an initial event of 'infinite' size, an event 
that spans the whole system. This is a reasonable situation because we use 
a uniform threshold distribution. In this case there is a finite probability to 
generate a stress 77 which exceedes the largest threshold, thus causing the whole 
system to reorganize. For some of the examples presented in Section 4 the 
probability to find an infinite event is even higher than 10 -5 . This probability 
can be considered relatively large in a system where one has to do about 
10 6 — 10 9 time-steps to get a good statistics. 



3. 1 Mean-field solution 

The exact way to calculate P t {s > S\) is the following. One has to compute 
the distribution pjif 2 '"'' T? *" 1 (a;) which is the distribution that arises if after the 
big event at time t = a sequence of stress values rji, r) 2 , ■ ■ ■ , r\ t _x is generated 
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during the following time steps. Then the equation 

xt(rii,V2,-,Vt-i,si) 

J pf™-^-\ x ')dx' = Sl (3) 



has to be solved. That gives the quantity x t (r)i,r)2, ■ ■ ■ , ?7t-i, Si), the thresh- 
old that has to be exceeded by the stress at time t to generate an event 
> s±. From the stress distribution one can then calculate the corresponding 
probability pj 11 ' 712 '---^- 1 ( s > Sl ). Finally the average over all possible sequences 
771 , 772, • • • , rjt-i has to be taken to end up with the exact solution for P t (s > si). 
Obviously there is no hope doing this analytically. 

Instead of the exact solution for P t (s > si) we can calculate a mean-field 
solution if we average over all possible sequences rji, i]2, ... , r) t - 1 before we solve 
Eq. (3). Note that in this context, the notion mean-field does not stand for the 
average state of the system, which does not tell us anything about aftershocks, 
but for the average fluctuations typically found in a time-intervall At. These 
average fluctuations are a measure for the return to the average state, after a 
large event has caused a significant departure from it. Consequently, the mean- 
field solution is time- dependent. For At — > 00, the time-dependent mean- field 
threshold distribution converges to the average threshold distribution, denoted 
by p(x) in [28] 

In Appendix A, we show that the averaging over all fluctuations in a time inter- 
vall of t time-steps equals to t times iterating the master-equation. Therefore, 
to calculate the mean-field solution for P t (s > si) we have to insert pt(x), the 
t-th iterate of the master-equation, into Eq. (3). The details of this calculation 
are given in Appendix B. 



3.2 Approximation for t 

In this paragraph we will calculate the dependency of the exponent r on s\ 
under the assumption that the probability to find aftershocks decays indeed 
as a power-law, i.e. that we can assume P t (s > si) ~ t~ T . A fairly simple 
argument shows that for the probability P t (s > si) to decrease as a power- law 
the exponent r must be proportional to 1 — si for s\ not too small. Again we 
concentrate on exponentially distributed stress only. 

We begin with an approximation of the quantity x t (si), which is the average 
threshold at time t above which a stress value must lie to trigger an event of 



6 



size > si. In the mean-field approximation, x t (si) is defined by the equation 

xt(si) 

J p t (x)dx = s 1 . (4) 
o 

Because pt{x) and s\ are normalized, we can rewrite this equation (again we 
assume pthresh(^) to be uniform between and 1): 

j p,(x)dx=l-s 1 . (5, 

xt(si) 

For the most reasonable stress distributions the distribution of the agents pt(x) 
forms a plateau in the region close to x — 1. Therefore for sufficient large s\ 
we can approximate the integral in Eq. (5) by substituting pt(x) with its value 
at x — 1, which is p<(l). Eq. (5) then becomes 

(l-x t (s 1 j)p t (l) = l-s 1 . (6) 

The values p t (l) are a function of t. We define 

R(t):= Pt (l) (7) 

and find for x t (si): 

si - 1 + R(t) 
*«(*) = — • W 

The probability P t (s > Si) now becomes 

P t (« > «i) = exp ( - — - j = exp ( - j . (9) 



The principal idea to derive a relation between r and si is as follows. The 
function R(t) is obviously independent of r and si. We make the ansatz P t (s > 
Si) ~ t~ r , rearrange Eq. (9) for R(t) and then get a condition on r and Si 
since they should cancel exactly. Hence we have to solve the equation 



at' 



where a is the constant of proportionality. We take the logarithm on both 
sides to get 

si - 1 + R(t) , s 
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and finally 



R(t) = — . (12) 

w 1 + alna -to hit v 1 



This is of the form Ci/(c 2 — lnt), where 



c = ^ (13) 



and 



1 + oTna 

c 2 = . (14) 

ra 

For every function of the form c\j (c 2 — lnt), the constants C\ and c 2 are unique, 
as can be seen if we write 

c 2 -lnt In sap' 1 ; 

A change in c 2 results in a rescaling of the variable t, while a change in c\ 
results in a rescaling of the whole function. Consequently, neither c\ nor c 2 can 
depend on r or s\. This can be seen as follows. If, e.g., c\ depended on s±, then 
a change in si would rescale the function R(t). But this function is independent 
of si according to its definition (Eq. (7)). Hence C\ must be independent of 
S\ in itself. A similar argument holds for the variable c 2 . Therefore, S\ and r 
must cancel exactly in Eq. (13). This leads to the condition 

r = i— (16) 



Up to now we have not done any assumptions about the size of the first big 
event after which we are measuring the subsequent aftershocks. Therefore the 
proportionality r ~ (1 — s\) should hold in general, as long as S\ is not too 
small. If we additionally assume the inital event to have infinite size (sq = 1) 
we can easily calculate the constant a in Eq. (10). The meaning of this constant 
is the probability to get an event of size > S\ immediately after the initial big 
event, as can be seen by setting t — 1: 

Pi(s > si) = ar T = a. (17) 

For the case s — 1 the distribution of thresholds po(x) is uniform and thus 

a = exp( ) = exp( ) . (18) 

a a 

With Eqs. (9), (10), (13), and (18) we can write the probability P t (s > si) as 

Pt(s > si) = e- Sl/a t- (1 - Sl)/(acl) . (19) 
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In Section 4 we will find numerically that c\ — a 1 , and therefore r = 1 — s±. 
3.3 Limiting cases 

For two limiting cases we can deduce the behaviour of the exponent r regard- 
less of the stress distribution. We begin with the case so — 1, s± — > 1. From 
Eq. (16) we find that r — > as si — > 1 under the assumption of an exponential 
stress distribution. But this result is more general. For s± — 1 the probability 
> s i) reads simply 



and hence is constant in time. Consequently we have r = for any stress 
distribution. From continuity we have r — > as si — > 1. 

A similar argument holds when either s or si approaches 0. For so = 0, the 
probability P t (s > si) reduces to the mean probability to find an event of 
size s > si. Hence r = 0. For s\ = 0, the probability Pt(s > si) becomes 1, 
because an event of size at least zero can be found in every time step. Hence 
also in this case r = 0. From continuity we have again r — > as either s — > 
or si — > 0. 



4 Numerical results 

In the previous section we have focused on the behaviour of the system in the 
aftermath of an infinite event. This situation is not only analytically tractable, 
but it also makes it simpler to obtain good numerical results. If we want to 
measure the probability to find aftershocks following events exceeding some 
finite but large s , we have to wait a long time for every single measurement 
since the number of those large events vanishes with a power-law. This makes it 
hard to get a good statistics within a reasonable amount of computing time. 
If, on the other hand, we focus on the situation of an infinite initial event, 
we can simply initialize the agents with the uniform threshold distribution, 
take the measurement up to the time t we are interested in and repeat this 
procedure until the desired accuracy is reached. Unless stated otherwise, the 
results reported below have been obtained in this way, and with exponentially 
distributed stress. 

The t-ih iteration of the master-equation should give exactly the average dis- 
tribution of the agent's thresholds at time t. In Fig. 1 it can be seen that this 
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(20) 



i 
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is indeed the case. The points, which represent simulation results, lie exactly 
on the solid lines, which stem from the exact analytical calculation. 



The mean-field approximation for P t (s > s±) should be valid if the agent's dis- 
tribution at time t does not fluctuate too much about the average distribution 
pt{x). Since there are many more small events than big ones the fluctuations 
should occur primarily in the region of small x. Consequently we expect the 
mean-field approximation to be valid for large si, and to break down for 
small si. Fig. 2 shows that already for moderately large si the mean-field 
approximation captures the behaviour of P t (s > si), with a slight tendency 
to underestimate the results of the measurement. Note that the statistics is 
becoming worse with increasing si due to the rapidly decreasing probability 
to find any events of size > si for large s±. 

In Fig. 3 a measurement of the probability P t (s > si) is presented for a number 
of simulations with different values of the parameter /. As it can be seen, the 
parameter / does not affect the exponent of the power-law, but limits the 
region where scaling can be observed. Note the difference between the effect 
seen here and typical cut off effects in the theory of phase transitions. The 
quantity P t (s > si) is not a probability distribution, but a time dependent 
probability, which tends towards a constant for t — > oo. Therefore, we do not 
see an exponential decrease at the cut off timescale. Instead, the probability 
Pt(s > si) levels out to the time- averaged value P(s > si), which is the 
average probability to find events of size s > s±. 

In section 3.2 we showed that r ~ 1 — si, under the condition of a sufficiently 
large s 1 . Simulations indicate that the constant of proportionality is just 1, 
which means the constant C\ in Eq. (13) equals a^ 1 . If this hypothesis is true, 
Eq. (19) becomes 



should yield a functional dependency P sca ied(^) = t -1 . Fig. 4 shows the results 
of such a rescaling for different a and s±. All the data-points lie exactly on 
top of each other in the region where the statistics is good enough (about 
t < 100). We find that for a up to 0.1, Eq. (21) is very accurate for si between 
about 0.1 and I. With smaller a, Eq. (21) holds even for much smaller si. 

The situation becomes more complicated if we study the sequence of after- 
shocks caused by an initial event of finite size. The probability to find an 
event of size s > si after some initial event of size s > s decreases with a 



Pt(s > Si) 



g-siA^-C 1 -*!) 



(21) 



This means, a rescaling of the form 




(22) 
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power-law, but the exponent is not a simple function of s\. Rather, it depends 
on So as well. In Fig. 6 we have displayed the results of a measurement with 
Si = 3 x 10~ 4 and several different Sq, ranging from 5 x 10~ 4 to 1. The curve for 
s = 1 has been obtained with the method described at the beginning of this 
section. We find that the aftershocks' decay pattern for s < 1 continuously 
approaches the one for so = 1 as so — > 1. This shows that it is indeed justified 
to study the system in the aftermath of an infinite initial event and then to 
extrapolate to finite but large initial events. Note that in Fig. 6, Si is so small 
that Eq. (21) does not hold anymore. 

Sneppen and Newman have argued that the decay pattern of the aftershocks 
is independent of the respective stress distribution. Our results do not support 
that. Despite the fact that the exponent of the power-law seems to be inde- 
pendent of a in the case of exponential stress, as we could show above, the 
exponent is not independent of the functional dependency of the stress distri- 
bution. If we impose, for example, gaussian stress with mean and variance 
a, we find (Fig. 5) exponents larger than 1 for moderate si. We do event find 
a qualitatively new behaviour of the system. The exponent is getting larger 
with increasing si, as opposed to the exponent getting smaller for exponential 
stress. Of course, this can only be true for intermediate s±. Ultimately, we 
must have r — > for s± — > 1. 

Finally, we present some results about systems with finite growth. In these sys- 
tems, there exists some competitive dynamics between the removal of agents 
with small thresholds through stress and their regrowth. Aftershocks appear 
in the infinite growth model because the reorganization event leaves a larger 
proportion of agents in the region of small thresholds, thus increasing the prob- 
ability for succeding large events. In the finite growth version, on the contrary, 
this can happen only if the regrowth is faster than the constant removal of 
agents with small thresholds through average stress. If the regrowth is too 
slow, the probability to find large events actually is decreased in the after- 
math of an initial large event. The interplay between these two mechanisms 
is shown in Fig. 7. The regrowth of species is done according to Eq. (2). For a 
small growth-rate g, the probability to find aftershocks is reduced significantly, 
and it aproaches the equilibrium value after about 100 time-steps. With larger 
g, the probability Pt(s > si) increases in time until a maximum well above the 
equilibrium level is reached, and then decreases again. The maximum moves 
to the left to earlier times t with increasing g. When g is so large that the 
maximum coincides with the post initial event, the original power-law is re- 
stored. Note that, as in the case of infinite growth, the measurement depends 
on the choice of the parameters so and s±. Consequently, with a different set 
of parameters the curves will look different, and the maximum will appear at 
a different time. Nevertheless, we find that the qualitative behaviour is not 
altered if we change s or s x . 
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Instead of logistic growth, we can also think of linear growth, i.e. AN = giV max , 
where g is again the growth rate. In order to keep the system finite, we stop 
the regrowth whenever the system size N exceedes the maximal system size 
iV max . In such a system, aftershocks can be seen for much smaller growth 
rates (Fig. 8). Note that apart from the growth rate, all other settings are 
identical in Fig. 7 and Fig. 8. Linear growth refills the system much quicker 
than logistic growth with the same growth rate. Therefore the time intervall 
in which aftershocks are suppressed is much shorter for linear growth. 



5 Conclusion 

We could show in the present paper that the decay pattern of the aftershock 
events depends on the details of the measurement. Although the qualitative 
features remain the same for different parameters s and si (e.g. a power- 
law decrease in the infinite-growth version), the quantitative features vary to 
a large extend. The exponent r of the power-law is significantly affected by 
an alteration of so or s±. Therefore the measurement-procedure proposed by 
Sneppen and Newman can reveal the complex structure of aftershock-events 
only if the change of the measured decay pattern with varying s and si is 
recorded over a reasonably large intervall of different values. This should be 
considered in a possible comparison between the aftershocks' decay pattern 
from a model and from experimental data. A more in-depth analysis could 
probably be achieved with the formalism of multifractality (see e.g. [33]). 

We found the aftershocks' decay pattern to vary with different stress distri- 
butions. This is in clear contrast to Sneppen and Newman. They reported a 
power-law with exponent t ~ 1 for the infinite growth version, independent of 
the respective stress distribution they used. The question remains why Snep- 
pen and Newman measured such an exponent in all their simulations. The 
answer to this lies in the fact that they did only simulations with so < 1. 
For the reasons explained at the beginning of Section 4, under this condition 
one has to choose a relatively small s , and accordingly, a very small s±. This 
causes the measurement almost inevitably to lie in the intermediate region 
between the limiting cases of Sec. 3.3. In this region, for the most reasonable 
stress distributions and a large array of different values for s and s±, we find 
indeed exponents around 1. 

The application of coherent-noise models to earthquakes has been discussed 
in [27]. Two very important observations regarding earthquakes, the Gutenberg- 
Richter law [34,35] and Omori's law [36], can be explained easily with a 
coherent-noise model. The Gutenberg-Richter law states that earthquake mag- 
nitudes are distributed according to a power-law. Omori's law, which interests 
us here, is a similar statement for the aftershocks' decay pattern of earth- 
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quakes. In the data from real earthquakes, the number of events larger than a 
certain magnitude Mi decreases as t~ T in the aftermath of a large initial earth- 
quake. Consequently, we can only apply infinite-growth coherent-noise models 
to earthquakes. But this is certainly no drawback, since we expect the thresh- 
olds against movement at various points along a fault (with which we identify 
the agents of the coherent-noise model) to reorganize almost instantaneously 
during an earthquake. 

The exponent r is not universal, but can vary significantly, e.g., in [1] from 
r = 0.80 to t = 1.58. This would cause problems if the statement was true 
that for coherent-noise models we have r ~ 1. But as we could show above, 
the exponent can assume values in the observed range, depending on the stress 
distribution, the size of the initial event, and the lower cut-off (which we called 
Mi for earthquakes and si throughout the rest of the paper). 

For a further comparison, it should be interesting to study the dependency of 
the exponent r on variation of the cut-off Mi in real data. We are only aware 
of a single work where that has been done [37]. Interestingly, the authors do 
not find a clear dependency t(M 1 ). Nevertheless, this is not a strong evidence 
against coherent-noise models, since the aftershock series analysed in [37] con- 
sists mainly of very large earthquakes with magnitude between 6 and about 
8, which does not allow a sufficient variation of M ± . Statistical variations in 
the exponent r are probably larger for this aftershock series than the possible 
variations because of an assumed t(Mi) dependency. 

Numerical simulations of the finite-growth version have revealed a much more 
complex structure of aftershock events than present in the infinite-growth ver- 
sion. The competition between regrowth of agents and agent removal through 
reorganization events leads to a pattern where the probability to find events 
after an initial large event is suppressed for short times, enhanced for inter- 
mediate times and then falls off to the background level for long times. This 
observation can be important for the application of coherent-noise models to 
biological extinction. It might be possible to identify a time of reduced and 
a time of enhanced extinction activity in the aftermath of a mass extinction 
event in the fossil record. This would be a good indication for biological ex- 
tinction to be dominated by external influences (coherent-noise point of view) 
rather than by coevolution (SOC point of view). 
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A Rederivation of the master-equation 



In this appendix we are interested in the average state a coherent-noise sys- 
tem will be found in several time-steps after some initial state with threshold 
distribution pt {x). Our calculations will lead to a rederivation of the master- 
equation for coherent-noise systems. Although a master-equation has already 
been given for the infinite-growth version and has been generalized to the 
finite-growth version, these master-equations have not been derived in a strin- 
gent way, but just have been written down intuitively. Our calculation will 
confirm the main terms of the previously used equations, but we will find an 
additional correcting term that becomes important for large /. 

Consider the case of infinite growth. At time t the threshold-distribution may 
be pt {x). We construct the distribution Pt +i(x), which is the distribution that 
arises at time to + 1 if a stress rj is generated at time to- A stress rj will cause 
a proportion = J^dx p to (x) of the agents to move. We have to distinguish 
two regions. For x < rj, all agents are removed. Then they are redistributed 
according to s,,Pthresh (x) ■ A small fraction / of the agents is then mutated, 
which results in 

Pt +l( X ) = (1 - f)SriPthT C sh(x) + /pthrcsh(^) 5 X < 1] . (A.l) 

For x > rj, the redistribution of the agents gives Pt (x) + s^Pthrcsh^)- With the 
subsequent mutation we obtain: 

Pt +l( X ) = I 1 - f)\Pto( X ) + S v Pthresh(x) ) + /Pthrcsh(a;) ] X > T] . 

V 7 (A.2) 



We take the average over i] to get the distribution p to+ i(x) that will on average 
be found one time-step after pt (x): 





oo ri 

= J dr) p stiess (r})p thiesh (x) (1 - /) J p to (x')dx' + f 



x 

stress 

(r))p t0 (x)(l - f) 



— Pthrcsh J 



/ + (!-/) / dr] p stress (v) / Pt a {x')dx' 



+ p to (x)(l - f)(l - p moye (x)) . 



(A.3) 



Here, p mov e( a; ) is the probability for an agent with threshold x to get hit 
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by stress, viz. p m0 ve(x) = S^ 3 dx' p stvcss {x') . To proceed further we have to 
interchange the order of integration in the remaining double integral. Note 
that f£°dr) J^dx' = f£°dx' f^dr), and therefore 



oo oo 



Pt +i(x) =Pthrcsh(a;) 



/ + (!-/) jdx' J dr}p stress (r})p to (x') 



+ p to (x)(l - f)(l - p moye (x)) 



Pthrcsh(a;) 



/ + (1 - /) J dx' p to (x')p movc (x') 



+ Pto(x)(l ~ /)(1 -Pmove(z)) 

oo 

Pthresh(^) J dx' (f + (1 - /)p I „ovc(a; / ))pto( a;/ ) 


+ p to (x)(l-/)(l-p move (x)) (A.4) 



We are thus led to the master-equation 

&Pt(x) = ( - / - p movc (x) + fpmovc(x) jPt(x) + Ap thrcsh (x) , 

V 7 (A.5) 

where A is the normalization constant f£°dx' (/ + (1 — f)p m ove(x')) p t (x') . 

We notice the appearance of the term fp movc (x) which was not present in the 
master-equation used by Sneppen and Newman. The term arises if one takes 
into account the fact that the agents which are hit by stress get new thresholds 
before the mutation takes place. Therefore every agent with threshold x has 
the probability fp m avc{x) to get two new thresholds in one time-step. But obvi- 
ously this is exactly the same as geting only one new threshold. Consequently, 
the term /p mov c( a; ) has to be present to avoid double-counting of those agents 
which are hit both by stress and by mutation. Nevertheless, this term does 
not affect the scaling behaviour of the system, because the derivation of the 
event size distribution in [28] has been done under the assumption /<1. 

Eq. (A.4) gives the average state of the system one time-step after the initial 
state pt (x). If we are interested in the average state t time-steps after the 
initial state, we have to repeat the calculations in Eqs. (A.1)-(A.4) t times. 
Since all averages in these calculations can be taken independently, this is 
exactly the same as t times iterating the master-equation (A.5). 
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B Calculation of the mean-field solution. 

We assume that at time t = a big event takes place and produces the distri- 
bution po(x). If we apply the master-equation (A. 5) t times to this distribution 
Po{x), we will end up with a distribution pt(x) that tells us the average state 
of the system at time t after the big event. 

In the following we will use 

T(x) :=(l-/)(l-?We(z)) (B.l) 

and write A t for the normalization constant that appears on the right-hand 
side of Eq. (A. 5) at time t. The average distribution at time t then becomes 

p t (x) = T(x)p t -i(x) + A t p thiesh (x) 

= T t (x) Po (x) + T t ' k (x)A kPthrcsh (x) . (B.2) 

k=l 

We integrate both sides of Eq. (B.2) and find a recursion relation for the 
constants A t : 

OO £ OO 

A t = 1 - f T t (x)p (x)dx + J2 A k ! T t - k (x)p thicsh (x)dx. 

o k=1 o (B.3) 

All integrals can be calculated analytically for a special choice of the thresh- 
old and stress distributions. As threshold distribution, we choose the uniform 
distribution Pthresh^) = 1; < x < 1, and as stress distribution we choose 
the exponential distribution p stre ss( ? 7) = exp(— r)/a)/a. Furthermore, we as- 
sume that the initial event was so large as to span the whole system, i.e. 
Po(x) = 1; < x < 1. 

Under the above assumptions there is basically one integral that appears in 
Eq. (B.3), which is 

i 

I n ■= J T n (x)dx 

o 
i 

= /(l-/r(l-e-^)V (B.4) 
o 

and Eq. (B.3) becomes 

t-i 

A t = l-I t + Y / h-kA k . (B.5) 

k=l 
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With the aid of the binomial expansion (1 + a) 11 = J2k=o u)a fc we find 




We are now in the position to calculate the probability that an event of size 
s > s\ occurs at time t after the initial big event. The minimal stress value 
7] m in that suffices on average to generate such an event is the solution to the 
equation 



Vm'm 



J dxp t (x) = s 1 . (B.7) 
o 

The corresponding probability is P t (s > Si) = exp(—r] min /a). We invert this 
expression and insert it into Eq. (B.7). The resulting equation determines the 
probability P t (s > Si): 

—a In Pt(s>si) 

J dxp t (x)= Sl . (B.8) 
o 

The integrals that appear after inserting pt{x) into the above equation are 
very similar to the integral /„ defined in Eq. (B.4). We define 

-crlnP 

J n (P):= J T n (x)dx. (B.9) 
o 

This integral can be taken in the same fashion as the calculation of /„ in 
Eq. (B.6). We find 



UP) = a(l -/)»(- In P + ± Q 1^(1 - P) 



(B.10) 



Eq. (B.8) now becomes 

Jt(Pt{s > sSj + J2 Jt-k(Pt(s > Sl )^A k = Sl . (B.ll) 

All the quantities which appear in this equation are given above in analytic 
form. Therefore solving Eq. (B.ll) is simply a problem of root-finding. With a 
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computer-algebra program such as Mathematica, the recursion relation for the 
constants Aj, as well as the sums that appear in the quantities I n and J n {P) 
can be evaluated analytically if we restrict ourselves to moderate t. Then the 
only numerical calculation involved in the computation of P t (s > si) is the 
calculation of the root of Eq. (B.ll). 
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Fig. 1. The average distribution of agents at time t = 1, t = 10, and t = 100 after the 
initial event of size oo. The solid line is the analytical result from the iteration of the 
master-equation, the dots show the simulation results. Parameters where a = 0.05 
and / = 10~ 5 
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Fig. 3. The probability Pt(s > s±) in simulations with a = 0.1, s\ = 0.06 and 
several different values for / (from bottom to top: / = 10~ 4 , / = 10~ 3 , / = 10~ 2 , 

/ = io- 1 ). 
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time t 

Fig. 4. The rescaled probability Pt(s > s\). The solid line shows a t~ x decrease for 
comparison. 
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Fig. 5. The probability Pt(s > s\) after an initial infinite event in a simulation 
with gaussian stress. Parameters are a = 0.1, / = 10~ 5 , and, from bottom to top, 
s\ = 0.2, si = 0.1, s\ = 0.02. In a simulation with gaussian stress the distribution 
is getting steeper with increasing si. 
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Fig. 6. The probability Pt(s > s±) after an initial finite event. Parameters are 
a = 0.05, / = 10~ 5 , si = 3 x 10~ 4 , and, from bottom to top, s = 5 x 10~ 4 , 
s = 2 x 10 -3 , s = 0.01, s = 0.1 and s = 1. 
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0.0001 1 ' — ' — 1 — ' — 

1 10 100 1000 10000 

time t 

Fig. 7. The probability Pt(si > 3 x 1CT 4 ) after an event of size so > 0.001 in a 
simulation with finite logistic growth. Parameters are a = 0.5, / = 10~ 5 , and, from 
bottom to top, g = 2 x 10~ 4 , g = 2 x 10~ 3 , g = 2 x 10~ 2 g = 2 x 10 _1 , g = 10. A 
power-law can only be seen for relatively large growth rates. For small growth-rates, 
the probability to find aftershocks is reduced significantly. 
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Fig. 8. The probability Pt(si > 3 x 10~ 4 ) after an event of size s > 0.001 in a 
simulation with finite linear growth. Parameters are a = 0.5, / = 10~ 5 , and, from 
bottom to top, g = 10~ 5 , g = xl0~ 4 , g = xl0~ 3 , g = 10. Aftershocks are seen for 
much smaller growth rates than in the version with logistic growth. 



23 



